In [1]:
from __future__ import division
from functions import *
from utils import *
from gevi_classes import *

%matplotlib inline
%load_ext autoreload
%autoreload 2
In [2]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 1000;
In [3]:
# instantiate utility class
gr = Graph()

Load DATA

Data collected in 5 X 1min, with 1min pause in between 2921 images obtained in 60s

In [4]:
username = os.path.expanduser('~').split('/')[-1]
if username == "GP1514":
    print("At Imperial")
    mouseAPath = '/Volumes/DATA/GEVI/MouseA/'
    mouseBPath = '/Volumes/DATA/GEVI/MouseB/'
    mouseM1217FPath = '/Volumes/DATA/GEVI/M1217F/'
    mouseM1223MPath = '/Volumes/DATA/GEVI/M1223M/'
else:
    print("Using laptop")
    mouseAPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 27/'
    mouseBPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 22/'
    mouseCPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 28/'

At Imperial
In [5]:
discard = loadDic()
In [6]:
mouseA = Mouse('mouseA', mouseAPath, [3,4,5,6],discard['MouseA'] )
# mouseA = Mouse('mouseA', mouseAPath, [3],discard['MouseA'] )
mouseB = Mouse('mouseB', mouseBPath, [2,3,4,5],discard['MouseB'] )

mouseM1217F = Mouse('mouseM1217F', mouseM1217FPath, range(11,14),discard['MouseC'] )
mouseM1223M = Mouse('mouseM1223M', mouseM1223MPath, range(35,36),discard['MouseC'] )
In [7]:
mouseA.loadData(SVD=False)
mouseB.loadData(SVD=False)
# mouseC.loadData()
Out[7]:
1
In [8]:
mouseB.experiments[1].repeats[0].setInfo('End of transition, then discard')
mouseB.experiments[2].repeats[0].setInfo('Transition to anesthesia')
mouseA.experiments[0].repeats[2].setInfo('Keep')
mouseA.experiments[2].repeats[1].setInfo('Transition to desynchronization: short periods of silence trigger recovery of hemodynamics')
mouseA.experiments[3].repeats[1].setInfo('Woke up - Discard')
mouseA.experiments[3].repeats[2].setInfo('Discard')
mouseA.experiments[3].repeats[3].setInfo('Discard')
mouseA.experiments[3].repeats[4].setInfo('Discard')
mouseA.getInfo()
mouseB.getInfo()

 - experiment 3 has 3 valid repeats: 
	 - Repeat 1 (not valid)
	 - Repeat 2 (valid)
	 - Repeat 3 (not valid): Keep
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
 - experiment 4 has 5 valid repeats: 
	 - Repeat 1 (valid)
	 - Repeat 2 (valid)
	 - Repeat 3 (valid)
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
 - experiment 5 has 5 valid repeats: 
	 - Repeat 1 (valid)
	 - Repeat 2 (valid): Transition to desynchronization: short periods of silence trigger recovery of hemodynamics
	 - Repeat 3 (valid)
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
 - experiment 6 has 4 valid repeats: 
	 - Repeat 1 (not valid)
	 - Repeat 2 (valid): Woke up - Discard
	 - Repeat 3 (valid): Discard
	 - Repeat 4 (valid): Discard
	 - Repeat 5 (valid): Discard
 - experiment 2 has 5 valid repeats: 
	 - Repeat 1 (valid)
	 - Repeat 2 (valid)
	 - Repeat 3 (valid)
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
 - experiment 3 has 3 valid repeats: 
	 - Repeat 1 (not valid): End of transition, then discard
	 - Repeat 2 (valid)
	 - Repeat 3 (valid)
	 - Repeat 4 (not valid)
	 - Repeat 5 (valid)
 - experiment 4 has 5 valid repeats: 
	 - Repeat 1 (valid): Transition to anesthesia
	 - Repeat 2 (valid)
	 - Repeat 3 (valid)
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
 - experiment 5 has 5 valid repeats: 
	 - Repeat 1 (valid)
	 - Repeat 2 (valid)
	 - Repeat 3 (valid)
	 - Repeat 4 (valid)
	 - Repeat 5 (valid)
In [9]:
# mouseA.experiments[0].repeats.mRatio

Plot Data

Mouse A

In [10]:
gr.plotHV(mouseA)

Mouse B

In [11]:
gr.plotHV(mouseB)

Mouse C

In [12]:
# gr.plotHV(mouseC)

Transfer functions

$V * \alpha = H$

$\mathscr{F}$ transform : $\hat{V}.\hat{\alpha} = \hat{H}$

$\alpha = \mathscr{F}^{-1} \frac{\hat{H}}{\hat{V}}$

ALPHA CURVE FITTING

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{\tau_3-\frac{t}{\tau_2}}$

Mouse A

In [13]:
def meanCorH(mouse):
    corH = np.mean(flatten(mouse.experiments.corH))
    return corH

def meanCorR(mouse):
    corR = np.mean(flatten(mouse.experiments.corR))
    return corR
    
# optimization for mouseA
def fn(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseB, alpha)
    corR = np.mean(flatten(mouseB.experiments.corR))
    res = (1 - corR)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def fnH(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseB, alpha)
    corH = np.mean(flatten(mouseB.experiments.corH))
    res = (1 - corH)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def computeAndPlotCorr(mouse, alpha):
        for i,exp in enumerate(mouse.experiments):
            for j,rep in enumerate(exp.repeats):
                gr.plotRModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.plotHModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.corrCoeff(mouse, exp = i, rep = j,alpha = alpha )
                
x = np.real(xax(mouseA.experiments[0].repeats[0].mRatio, 60000))
In [14]:
r= mouseA.experiments[0].repeats[0].ratio
out=np.ones((60,60))
In [35]:
# filtering parameters
mouseA.minFreqAlpha = 1
mouseA.maxFreqAlpha = 200
mouseB.minFreqAlpha = 1
mouseB.maxFreqAlpha = 200
mouseA.window = 10
gr.fHmax = 30
gr.fRmax = 200

Tranfer functions mouse A

In [36]:
gr.plotTF(mouseA)
Out[36]:
1
In [37]:
gr.plotmTF(mouseA)
Out[37]:
1

Transfer functions mouse B

In [38]:
gr.plotTF(mouseB)
Out[38]:
1
In [39]:
gr.plotmTF(mouseB)
Out[39]:
1

Voltage and Hemo reconstruction : with mean alpha function

In [40]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
alpha = mouseB.experiments[0].repeats[2].alpha
gr.plotRModel(mouseA, exp = 2, rep = 4,alpha = alpha )
gr.plotHModel(mouseA, exp = 2, rep = 4,alpha = alpha )

(-4.42534470505-547.440984603j)
In [41]:
alpha = mouseB.experiments[0].repeats[2].meanAlphaModel
gr.plotRModel(mouseA, exp = 2, rep = 4,alpha = alpha )
gr.plotHModel(mouseA, exp = 2, rep = 4,alpha = alpha )

(749.403239322-561.685929646j)

Voltage and Hemo reconstruction : with mean alpha model function

In [42]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
alpha = mouseB.experiments[2].meanAlphaModel
gr.plotRModel(mouseA, exp = 2, rep = 4,alpha = alpha )
gr.plotHModel(mouseA, exp = 2, rep = 4,alpha = alpha )

(-112.232725735-204.364733704j)

Distribution of alpha model parameters $\tau_0, \tau_1, \tau_2, \tau_3$

$v * \alpha = h$

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{-\frac{t}{\tau_2}}$

In [43]:
alpha = mouseB.experiments[2].meanAlphaModel
fontsize=16
plt.plot(xax(alpha,60000),alpha)
plt.xlabel('Time [s]')
plt.title('Alpha model of mean TF mouse B')
plt.text(30,-0.01, r'$ \tau_{0} = %.2g$' %mouseB.experiments[2].meanAlphaParams[0], fontsize =fontsize)
plt.text(30,-0.02, r'$ \tau_{1} = %.2g$' %mouseB.experiments[2].meanAlphaParams[1], fontsize =fontsize)
plt.text(30,-0.03, r'$ \tau_{2} = %.2g$' %mouseB.experiments[2].meanAlphaParams[2], fontsize =fontsize)
Out[43]:
<matplotlib.text.Text at 0x11ef617f0>
In [44]:
gr.computeCorr(mouseA, alpha=None)
gr.computeCorr(mouseB, alpha=None)
gr.computeAlphaModels(mouseA)
gr.computeAlphaModels(mouseB)
gr.plotParamsIndex(mouseA)
gr.plotParamsIndex(mouseB)
gr.plotParamsCor(mouseA)
gr.plotParamsCor(mouseB)
print([meanCorR(mouseA),meanCorH(mouseA),meanCorR(mouseB), meanCorH(mouseB)])

[0.23142305115075679, 0.27323145493387258, 0.26343748355642055, 0.24467684575888901]

Optimize for hemo

In [45]:
# ?minimize
In [46]:
resultHemo = minimize(fnH, x0=[1,1,1], method = 'COBYLA',
                      options={'gtol': 1e-8, 'disp': True})
print(resultHemo)

     fun: 0.40569510861206187
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 122
  status: 1
 success: True
       x: array([ 1.19414028,  0.54609966,  1.64826181])
In [47]:
p_est = resultHemo.x
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
plt.plot(xax(alpha,60000),alpha)
# alpha2 = mouseB.experiments[2].meanAlphaModel
# plt.plot(xax(alpha2,60000),alpha2)
Out[47]:
[<matplotlib.lines.Line2D at 0x11f423e10>]
In [48]:
# computeAndPlotCorr(mouseA, alpha)
gr.plotRModel(mouseA,2,4,alpha)
gr.plotHModel(mouseA,2,4,alpha)
# gr.plotHModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
# gr.plotRModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
def printl(tab):
    for t in tab:
        print(t)

(538.808174529-407.424321391j)

Optimize for voltage

In [49]:
# bounds = [(-10,10),(-10,10),(-10,10),(-10,10)]
# resultVolt = differential_evolution(fn, bounds, init = 'random')
resultVolt = minimize(fn, method = 'COBYLA', x0=[-1e-2,1,10,1],
                      options={'gtol': 1e-6, 'disp': True})
print(resultVolt)

     fun: 0.50934167175911094
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 438
  status: 1
 success: True
       x: array([ 0.4121683 ,  1.09384498,  3.68435018,  0.98029149])
In [50]:
p_est = resultVolt.x
print(p_est)
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]

[ 0.4121683   1.09384498  3.68435018  0.98029149]
In [51]:
plt.title('Optimized Alpha model of mean TF mouse B, exp %d'%(2+mouseB.start))
plt.text(30,0.5, r'$ \tau_{0} = %.2g$' %p_est[0], fontsize =fontsize)
plt.text(30,1.5, r'$ \tau_{1} = %.2g$' %p_est[1], fontsize =fontsize)
plt.text(30,2.5, r'$ \tau_{2} = %.2g$' %p_est[2], fontsize =fontsize)
plt.plot(xax(alpha,60000),alpha)
Out[51]:
[<matplotlib.lines.Line2D at 0x11d1636a0>]
In [52]:
gr.plotRModel(mouseA,2,4,alpha)
gr.plotHModel(mouseA,2,4,alpha)

(-127.26852031-649.22868781j)